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Abstract 

The Support Vector Machine (SVM) is a new and very promising classification technique developed by Vapnik and 
his group at AT&T Bell Laboratories [3, 6, 8, 24]. This new learning algorithm can be seen as an alternative training 
technique for Polynomial, Radial Basis Function and Multi-Layer Perceptron classifiers. The main idea behind the 
technique is to separate the classes with a surface that maximizes the margin between them. An interesting property 
of this approach is that it is an approximate implementation of the Structural Risk Minimization (SRM) induction 
principle [23]. The derivation of Support Vector Machines, its relationship with SRM, and its geometrical insight, 
are discussed in this paper. 

Since Structural Risk Minimization is an inductive principle that aims at minimizing a bound on the generalization 
error of a model, rather than minimizing the Mean Square Error over the data set (as Empirical Risk Minimization 
methods do), training a SVM to obtain the maximum margin classifier requires a different objective function. This 
objective function is then optimized by solving a large-scale quadratic programming problem with linear and box 
constraints. The problem is considered challenging, because the quadratic form is completely dense, so the memory 
needed to store the problem grows with the square of the number of data points. Therefore, training problems 
arising in some real applications with large data sets are impossible to load into memory, and cannot be solved 
using standard non-linear constrained optimization algorithms. 

We present a decomposition algorithm that can be used to train SVM's over large data sets. The main idea 
behind the decomposition is the iterative solution of sub-problems and the evaluation of, and also establish the 
stopping criteria for the algorithm. We present previous approaches, as well as results and important details of our 
implementation of the algorithm using a second-order variant of the Reduced Gradient Method as the solver of the 
sub-problems. 

As an application of SVM's, we present preliminary results in Frontal Human Face Detection in images. This 
application opens many interesting questions and future research opportunities, both in the context of faster and 
better optimization algorithms, and in the use of SVM's in other pattern classification, recognition, and detection 
applications. 
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1 Introduction 

In this report we address the problem of implementing a new pattern classification technique recently 
developed by Vladimir Vapnik and his team at AT&T Bell Laboratories [3, 6, 8, 24] and known as Support 
Vector Machine (SVM). SVM can be thought as an alternative training technique for Polynomial, Radial 
Basis Function and Multi-Layer Perceptron classifiers, in which the weights of the network are found by 
solving a Quadratic Programming (QP) problem with linear inequality and equality constraints, rather 
than by solving a non-convex, unconstrained minimization problem, as in standard neural network training 
techniques. Since the number of variables in the QP problem is equal to the number of data points, when 
the data set is large this optimization problem becomes very challenging, because the quadratic form is 
completely dense and the memory requirements grow with the square of the number of data points. We 
present a decomposition algorithm that guarantees global optimality, and can be used to train SVM's over 
very large data sets (say 50,000 data points). The main idea behind the decomposition is the iterative 
solution of sub-problems and the evaluation of optimality conditions which are used both to generate 
improved iterative values, and also establish the stopping criteria for the algorithm. 

We demonstrate the effectiveness of our approach applying SVM to the problem of detecting frontal faces 
in images, which involves the solution of a pattern classification problem (face versus non-face patterns) 
with a large data base (50,000 examples). The reason for the choice of the face detection problem as an 
application of SVM is twofold: 1) the problem has many important practical applications, and received a 
lot of attention in recent years; 2) the difficulty in the implementation of SVM arises only when the data 
base is large, say larger than 2,000, and this problems does involve a large data base. 

The paper is therefore divided in two main parts. In the first part, consisting of sections 2 and 3 we 
describe the SVM approach to pattern classification and our solution of the corresponding QP problem in 
the case of large data bases. In the second part (section 4) we describe a face detection system, in which 
the SVM is one of the main components. In particular, section 2 reviews the theory and derivation of 
SVM's, together with some extensions and geometrical interpretations. Section 3 starts by reviewing the 
work done by Vapnik et al. [5] in solving the training problem for the SVM. Section 3.1.1 gives a brief 
description and references of the initial approaches we took in order to solve this problem. Section 3.2 
contains the main result of this paper, since it presents the new approach that we have developed to solve 
Large Database Training problems of Support Vector Machines. 

Section 4 presents a Frontal Human Face Detection System that we have developed as an application of 
SVM's to computer vision object detection problems. 

2 Support Vector Machines 

In this section we introduce the SVM classification technique, and show how it leads to the formulation 
of a QP programming problem in a number of variables that is equal to the number of data points. We 
will start by reviewing the classical Empirical Risk Minimization approach, and showing how it naturally 
leads, through the theory of VC bounds, to the idea of Structural Risk Minimization (SRM), which is a 
better induction principle, and how SRM is implemented by SVM. 

2.1 Empirical Risk Minimization 

In the case of two-class pattern recognition, the task of learning from examples can be formulated in the 
following way: given a set of decision functions 

{/a(x):AgA}, /a: ^^{_i,i } 

where A is a set of abstract parameters, and a set of examples 

(xi,2/i),...,(x^,%), x; G s R N ,y l e {-1,1} 
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drawn from an unknown distribution P(x,y), we want to find a function f\* which provides the smallest 
possible value for the expected risk: 



RW = J\h{*)-y\P{*,y)dxdy 



The functions f\ are usually called hypothesis, and the set {/\( x ) : A G A} is called the hypothesis space 
and denoted by 7i. The expected risk is therefore a measure of how good an hypothesis is at predicting the 
correct label y for a point x. The set of functions f\ could be for example, a set of Radial Basis Functions 
or a Multi-Layer Perceptron with a certain number of hidden units. In this case, the set A is the set of 
weights of the network. 

Since the probability distribution P(x, y) is unknown, we are unable to compute, and therefore to minimize, 
the expected risk P(A). However, since we have access to a sampling of P(x, y), we can compute a stochastic 
approximation of P(A) , the so called empirical risk: 

1 l 

Pemp(A) = - ^ I/a(x 8 ) - Vi\ 
1 i=l 

Since the law of large numbers guarantees that the empirical risk converges in probability to the expected 
risk, a common approach consists in minimizing the empirical risk rather than the expected risk. The 
intuition underlying this approach (the Empirical Risk Minimization Principle) is that if P e mp converges 
to P, the minimum of P em p may converge to the minimum of P. If convergence of the minimum of P em p 
to the minimum of P does not hold, the Empirical Risk Minimization Principle does not allow us to make 
any inference based on the data set, and it is therefore said to be not consistent. As shown by Vapnik 
and Chervonenkis [25, 26, 23] consistency takes place if and only if convergence in probability of P em p 
to P is replaced by uniform convergence in probability. Vapnik and Chervonenkis [25, 26, 23] showed 
that necessary and sufficient condition for consistency of the Empirical Risk Minimization Principle is the 
finiteness of the VC-dimension h of the hypothesis space 7i. The VC-dimension of the hypothesis space 7i 
(or VC-dimension of the classifier f\) is a natural number, possibly infinite, which is, losely speaking, the 
largest number of data points that can be separated in all possible ways by that set of functions f\. The 
VC-dimension is a measure of the complexity of the set 7i, and it is often, but not necessarily, proportional 
to the number of free parameters of the classifier f\. 

The theory of uniform convergence in probability developed by Vapnik and Chervonenkis also provides 
bounds on the deviation of the empirical risk from the expected risk. A typical uniform Vapnik and 
Chervonenkis bound, which holds with probability 1 — tj, has the following form: 



P(A) < Pemp(A) 



\ 



Mmf + 1 -In | 

■^—^ VA G A (1) 



where h is the VC-dimension of f\. From this bound it is clear that, in order to achieve small expected risk, 
that is good generalization performances, both the empirical risk and the ratio between the VC-dimension 
and the number of data points has to be small. Since the empirical risk is usually a decreasing function 
of h, it turns out that, for a given number of data points, there is an optimal value of the VC-dimension. 
The choice of an appropriate value for h (which in most techniques is controlled by the number of free 
parameters of the model) is crucial in order to get good performances, especially when the number of data 
points is small. When using a Multilayer Perceptron or a Radial Basis Functions network, this is equivalent 
to the problem of finding the appropriate number of hidden units. This problem is known to be difficult, 
and it is usually solved by some sort of cross-validation technique. 

The bound (1) suggests that the Empirical Risk Minimization Principle can be replaced by a better 
induction principle, as we will see in the next section. 




2.2 Structural Risk Minimization 

The technique of Structural Risk Minimization developed by Vapnik [23] is an attempt to overcome the 
problem of choosing an appropriate VC-dimension. It is clear from eq. (1) that a small value of the 
empirical risk does not necessarily imply a small value of the expected risk. A different induction principle, 
called the Structural Risk Minimization Principle, has been proposed by Vapnik [23]. The principle is 
based on the observation that, in order to make the expected risk small, both sides in equation (1) should 
be small. Therefore, both the VC-dimension and the empirical risk should be minimized at the same time. 
In order to implement the SRM principle one needs then a nested structure of hypothesis spaces 

rdcn.2 c...cn n c... 

with the property that h(n) < h(n + 1) where h(n) is the VC-dimension of the set 7i n . Then equation (1) 
suggests that, disregarding logarithmic factors, the following problem should be solved: 

min R emp [\] 

rin \ 

The SRM principle is clearly well founded mathematically, but it can be difficult to implement for the 
following reasons: 

1. The VC-dimension of 7i n could be difficult to compute, and there are only a small number of models 
for which we know how to compute the VC-dimension. 

2. Even assuming that we can compute the VC-dimension of TC n , it is not easy to solve the minimization 
problem (2). In most cases one will have to minimize the empirical risk for every set TC n , and then 
choose the 7i n that minimizes eq. (2). 

Therefore the implementation of this principle is not easy, because it is not trivial to control the VC- 
dimension of a learning technique during the training phase. The SVM algorithm achieves this goal, 
minimizing a bound on the VC-dimension and the number of training errors at the same time. In the 
next section we discuss this technique in detail, and show how its implementation is related to quadratic 
programming. 

2.3 Support Vector Machines: Mathematical Derivation 

In this section we describe the mathematical derivation of the Support Vector Machine (SVM) developed 
by Vapnik [24]. The technique is introduced by steps: we first consider the simplest case, a linear classifier 
and a linearly separable problem; then a linear classifier and non-separable problem, and finally a non-linear 
classifier and non-separable problem, which is the most interesting and useful case. 

2.3.1 Linear Classifier and Linearly Separable Problem 

In this section we consider the case in which the data set is linearly separable, and we wish to find the 
"best" hyperplane that separates the data. For our purposes, linearly separable means that we can find a 
pair (w, b) such that: 



w • x 8 - + b > 1 Vx 8 - G Class 1 (3) 

w • x; + b < -1 Vx; G Class 2 (4) 

The hypothesis space in this case is therefore the set of functions given by 

fw,b = sign(w • x + b) (5) 
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Notice that if the parameters w and 6 are scaled by the same quantity, the decision surface given by (5) 
is unchanged. In order to remove this redundancy, and to make each decision surface correspond to one 
unique pair (w,6), the following constraint is imposed: 

min |w • x,- + 61 = 1 (6) 

i=i,...,e 

where xi, . . . ,x^ are the points in the data set. The set of hyperplanes that satisfy (6) are called Canonical 
Hyperplanes. Notice that all linear decision surfaces can be represented by Canonical Hyperplanes, and 
constraint (6) is just a normalization, which will prove to be very convenient in the following calculations. 
If no further constraints are imposed on the pair (w, 6) the VC-dimension of the Canonical Hyperplanes is 
N + 1 [24], that is, the total number of free parameters. In order to be able to apply the Structural Risk 
Minimization Principle we need to construct sets of hyperplanes of varying VC-dimension, and minimize 
both the empirical risk (the training classification error) and the VC-dimension at the same time. A 
structure on the set of canonical hyperplanes is defined by constraining the norm of the vector w. In fact, 
Vapnik shows that, if we assume that all the points xi, . . . , x^ lie in the unit N- dimensional sphere, the set 

{fw,b = sign(w • x + 6) | ||w|| < A} (7) 

has a VC-dimension h that satisfies the following bound [24] [23]: 

h < min{ \A 2 ] , N} + 1 (8) 

If the data points lie inside a sphere of radius R, then (8) becomes h < mm{\R 2 A 2 ~\, N} + 1. The 
geometrical reason for which bounding the norm of w constraints the set of canonical hyperplanes is very 
simple. It can be shown that the distance from a point x to the hyperplane associated to the pair (w,6) 
is: 

rf(x;w,6)=J — „ ' (9) 

||w|| 

According to the normalization (6) the distance between the canonical hyperplane (w,6) and the closest 
of the data points is simply tt^tt. Therefore, if ||w|| < A then the distance of the canonical hyperplane to 

the closest data point has to be larger than -j. We can then conclude that the constrained set of canonical 
hyperplanes of eq. (7) is the set of hyperplanes whose distance from the data points is at least -j. This is 
equivalent to placing spheres of radius -j around each data point, and consider only the hyperplanes that 
do not intersect any of the spheres, as shown in figure (1). 

If the set of examples is linearly separable, the goal of the S VM is to find, among the Canonical Hyperplanes 
that correctly classify the data, the one with minimum norm, or equivalently minimum ||w|| 2 , because 
keeping this norm small will also keep the VC-dimension small. It is interesting to see that minimizing 
||w|| 2 (in this case of linear separability) is equivalent to finding the separating hyperplane for which 
the distance between the two convex hulls (of the two classes of training data), measured along a line 
perpendicular to the hyperplane, is maximized. In the rest of this paper, this distance will be referred to 
as the margin. Figure (2) gives some geometrical interpretation of why better generalization is expected 
from a separating hyperplane with large margin. 

To construct the maximum margin or optimal separating hyperplane, we need to correctly classify the 
vectors x 8 - of the training set 

(xi,yi ),..., (*e,ye), x; e $l N 

into two different classes yi £ { — 1, 1}, using the smallest norm of coefficients. This can be formulated as 
follows: 

5 




Figure 1: Bounding the norm of w is equivalent to constraint the hyperplanes to remain outside spheres 
of radius -j centered around the data points. 



Minimize <I>(w) = — llwll 

V ) 2 n II 

w, b 

subject to 

yi(wxi + b) >1 i = l...t 

(10) 

At this point, this problem can be solved using standard Quadratic Programming (QP) optimization 
techniques and is not very complex since the dimensionality is JV + 1. Since N is the dimension of the 
input space, this problem is more or less tractable for real applications. Nevertheless, in order to easily 
explain the extension to nonlinear decision surfaces (which will be described in section 2.3.3), we look at 
the dual problem, and use the technique of Lagrange Multipliers. We construct the Lagrangian 

1 ' 



Z(w, 6, A) = -||w|| 2 - J2 Mw(w • x,- + b) - 1], (11) 



8 = 1 



where A = (Ai, . . . , A^) is the vector of non-negative Lagrange multipliers corresponding to the constraints 
in (10). 

The solution to this optimization problem is determined by a saddle point of this Lagrangian, which has 
to be minimized with respect to w and 6, and maximized with respect to A > 0. Differentiating (11) and 
setting the results equal to zero we obtain: 



dL(w,b,A) 



w 



<9w . , 

8 = 1 



X>«-W x «- = ° (12) 




Figure 2: (a) A Separating Hyperplane with small margin, (b) A Separating Hyperplane with larger 
margin. A better generalization capability is expected from (b). 



flZ(w,6,A) 
db 



Yl X ^ = ° 



:i3i 



8 = 1 



Using the superscript * to denote the optimal values of the cost function, from equation (12) we derive: 



w 



J2 A ^ x « 



14) 



8 = 1 



which shows that the optimal hyperplane solution can be written as a linear combination of the training 
vectors. Notice that only those training vectors x 8 - with A 8 - > contribute in the expansion (14). 
Substituting (14) and (13) into (11) we obtain: 



1 



£ £ 



F ( A ) = J2 X <- oll w l 2 = J2 A « - o J2 J2 A ' A jK» x i • x 



(15) 



8 = 1 



8=1 



8 = 1 j = l 



Writing (15) in matrix notation, incorporating non-negativity of A and constraint (13) we get the following 
dual quadratic program: 

Maximize F(A) = A • 1 - -A • DA 

subject to 

A-y =0 

A > 



where y = (yi, . . . , yi) and D is a symmetric t X t matrix with elements Dij = yiyj~x.i ■ Xj. 
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Notice that complementary slackness conditions of the form: 

A t %(w* • x t - + &*) - 1] = i = l,...,£ (16) 

imply that A 8 - > only when constraint (10) is active. The vectors for which A 8 - > are called Support 
Vectors. From equation (16) it follows that b* can be computed as: 

b = y % - w • x; 

for any support vector x 8 -. By linearity of the dot product and equation (14), the decision function (5) can 
then be written as: 

/(x) = sign \J2yiK(x-*i) + A (17) 

2.3.2 The Soft Margin Hyperplane: Linearly Non-Separable Case 

We now consider the case in which we still look for a linear superating surface, but a separating hyperplane 
does not exist, so that it is not possible to satisfy all the constraints in problem (10). In order to deal 
with this case one introduces a new set of variables {£i}f =1 , that measure the amount of violation of the 
constraints. Then the margin is maximized, paying a penalty proportional to the amount of constraint 
violations. Formally, one solves the following problem: 



18] 



i = l,...,£ (19) 

i = l,...,£ (20) 

where C and k are parameters which have to be determined beforehand and define the cost of constraints 
violation. Other monotonic convex functions of the errors can be defined (see [8] for the more general case). 
Notice that minimizing the first term in (18) amounts to minimizing the VC-dimension of the learning 
machine, thereby minimizing the second term in the bound (1). On the other hand, minimizing the second 
term in (18) controls the empirical risk, which is the first term in the bound (1). This approach, therefore, 
constitutes a practical implementation of Structural Risk Minimization on the given set of functions. In 
order to solve problem (18), we construct the Lagrangian: 

, i i i 

i(w,6, A,H,T) = -||w|| 2 - $>.-[2ft(w • x t - + b) - 1 + &] -$>.■& + C(£ii)\ (21) 

i=l i=l i=l 

where the non-negative multipliers A = (Ai, . . . , A^) and T = (71, . . . , 7^) are associated with constraints 
(19) and (20) respectively. The solution to this problem is determined by the saddle point of this La- 
grangian, which has to be minimized with respect to w, H and 6, and maximized with respect to A > 
and r > 0. Differentiating (21) and setting the results equal to zero, we obtain: 

0Z(w,6,A,E,r) ^ 

= (w- 2^Ai2/iXi) = (22) 



Minimize $(w, H) = — | 


|w|| 2 + C(f>) fc 

2 = 1 


w,6,H 




subject to 




y t (w -Xi + b) > 1 - 


-& 


& >o 





M(W '^ g - r| = X>,,, = » ,23, 



db 



ax(w,&,A,a,r) = I fcc (Eii e.j - a 8 - 7i = o * > 1 

dS 1 C-A,-- 7t - = k = l. 



8 = 1 
fc-1 



When k > 1, by denoting 



vCA; 

8 = 1 



we can rewrite equation (24) as: 



From equation (22) we obtain: 



124) 



x> = i^r. <*> 



tf - A,- - 7 8 = 0. (26) 



E A ^' x « ( 27 ) 



w = 

8=1 

Substituting (27), (23) and (25) into (21) we obtain: 

F(A, «) = E^'-EE WjViyjXi ■ x, ^- 1 - -r (28) 

8=i 8=1 j=i (&C)*-i v fc/ 

Therefore, in order to obtain the Soft Margin separating hyperplane we solve: 

k 

Maximize F(A, 6) = A • 1 - ±A • DA ^- ( 1 - I 

1 (kc)— v k 

subject to 

A-y =o 

A <S1 

A > 



(29) 



where y = (yi, . . . , y$) and D is a symmetric t X t matrix with elements Dij = yiyj~x.i ■ Xj. 

When k = 1, that is, penalizing linearly the violations in constraint (19), the set of equations (29) simplifies 

to: 

Maximize F(A) = A ■ 1 - \A ■ DA 
subject to 

A-y =0 (30) 

A <C1 

A > 

The value k = 1 is assumed for the rest of this paper, since it simplifies the mathematical formulation and 
has shown very good results in practical applications. By the linearity of the dot product and equation 
(27), the decision function (5) can be written as: 



/(x) = sign Y^ %'A*(x • x; 

\8=1 



b* (31) 



where b* = yi — w* • x 8 -, for any support vector x 8 - such that < \ < C (that is a support vector which is 
correctly classified). In order to verify this, notice that complementary slackness in the conditions of the 
form: 

X*[y t (w* • x,- + b*) - 1 + £] = i = l,...,£ (32) 

imply that A 8 - > only when constraint (19) is active, establishing the need for A 8 - > 0. On the other hand, 
(19) can be active due to £,- > 0, which is not acceptable since x 8 - would be a misclassified point. For k = 1 
in equation (24) we have 7; = C — A,-. Since 7; is the multiplier associated with constraint (20), 7; > 
implies £,- = 0, establishing the sufficiency of A 8 - < C. Notice that this is a sufficient condition, since both 
7i and A 8 - could be equal to zero. 

Note: 

Our calculation above of the threshold value b assumes the existence of some A 8 - such that < A 8 - < C. 
We have not found a proof yet of the existence of such A,-, or conditions under which it does not exist. 
However, we think this is a very reasonable assumption, because it is equivalent to the assumption that 
there is at least one support vector which is correctly classified. So far our computational results indicate 
that this assumption is correct, and we will use it in the rest of this paper. 

2.3.3 Nonlinear Decision Surfaces 

Previous sections have only treated linear decision surfaces, which are definitely not appropriate for many 
tasks. The extension to more complex decision surfaces is conceptually quite simple, and is done by mapping 
the input variable x in a higher dimensional feature space, and by working with linear classification in that 
space. More precisely, one maps the input variable x into a (possibly infinite) vector of "feature" variables: 

x -► <£(x) = (ai^(x), a 2 ^ 2 (x), . . . , a n <f> n (x), . . .) (33) 

where {a n }^ =1 are some real numbers and {^> n }^ 1 are some real functions 1 . The Soft Margin version of 
SVM is then applied, substituting the variable x with the new "feature vector" </>(x). Under the mapping 
(33) the solution of a SVM has the form: 



/(x) = sign (0(x) • w* + b*) => sign ^ wA,>(x) • #(x,- 



b* (34) 



A key property of the SV machinery is that the only quantities that one needs to compute are scalar 
products, of the form </>(x) • <j>(y)- It is therefore convenient to introduce the so-called kernel function K: 



A'(x, y) = <Kx) • cf>(y) = J2 a? n M*)My) (35) 

71 = 1 

Using this quantity the solution of a SVM has the form: 

/(x) = sign \J2 Vi A* # (x, x,-) + b* J (36) 

and the quadratic programming problem (30) becomes: 



The numbers {a„}™ =1 are clearly unnecessary, and could be absorbed in the definition of the {4> n }%Li, but we use them 
here just because they make the formulation easier. 
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Maximize F(A) = A • 1 - ±A • DA 
subject to 

A-y =0 (37) 

A <C1 

A > 

where D is a symmetric, semi-positive definite, £ X £ matrix with elements Dij = y 8 yjli'(x 8 , Xj). Notice 
that the decision surface (36) is now a nonlinear function, given by linear superposition of kernel functions, 
one for each support vector. The idea of expanding the input space in a feature space is therefore useful 
only if we find some solution to the following problem: starting from the feature space or starting from the 
kernel. 

Problem 2.1 Find a set of coefficients {a, n }^L 1 and a set of features {<j), n }^L 1 such that: 

1. the scalar product ii'(x, y) = <?!>(x) • <?!>(y) is well defined (for example the series converges uniformly); 

2. the scalar product ii'(x, y) = <?!>(x) • <?!>(y) is easy to compute as a function of x and y; 

In addition to these requirements, we also should require the features (f>j to be such that the scalar product 
ii'(x, y) defines a class of decision surfaces which is "rich" enough (for example includes some well-known 
approximation schemes). There are two different approaches to this problem. 

Starting from the feature space 

One approach consists in choosing carefully a set of features with "good" properties. For example, an 
obvious choice would be to take as features <^ 8 (x) monomials in the variable x up to a certain degree. 
Assuming, for simplicity, to work in a one- dimensional space, one could choose: 

(f)(x) = (1, x, x 2 , . . . , X ) 

where d could be very high, and the coefficients a 8 - are all equal to one. In this case the decision surface 
is linear in the components of <j>, and therefore a polynomial of degree d in x. This choice is unfortunate, 
however, because the scalar product 

<f>(x) ■ <f>(y) = l + xy + (xyf + ... (xyf 

is not particularly simple to compute when d is high. However, it is easy to see that, with a careful choice 
of the parameters a 8 - things simplify. In fact, choosing 

an = ( i ) 



it is easy to see that 



<t>(x).<t>(y) = Y,{ 

n=0 \ 



d n {xvY = (i + *y) d 



which considerably reduces the computation. A similar result, although with a more complex structure of 
the coefficients a n , is true in the multivariable case, where the dimensionality of the feature space grows 
very quickly with the number of variables. For example, in two variables we can define: 

<£(x) = (1, V% x 1 ,V2 x 2 ,xl,x 2 ,,V2 x 1 x 2 ) (38) 

In this case it is easy to see that: 
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ii-(x,y) = <£(x)-4>(y) = (l + x-y) 2 (39) 

It is straightforward to extend this example to the d- dimensional case. For example, in 3 dimensions we 
have: 

<£(x) = (l,V^Il,>/2l2,V^l3,Il,S2! :I; 3!V^Il2;2,V^Ill3,V2l2S3) 

and the scalar product is still of the form of eq. (39). Still in 2 dimension we can use features which are 
monomials of degree 3: 

</>(x) = (1, v3 x\, v3 x 2 , v3 xf, v3 x 2 ,, vGxix 2 , v3 xfx 2 , v3 x\x\, x\ , x\ ) 



and it can be shown that: 



ii'(x,y)=(l + x-y) 



It can also be shown that if the features are monomials of degree less or equal to d, it is always possible to 
find numbers a n in such a way that the scalar product is 

ii'(x,y)=(l + x-y) d (40) 

In the following we provide a few more examples of how one could choose the features first, and then, with 
a careful choice of the coefficients a n , arrive at an analytical expression for the kernel K. 

Infinite dimensional feature spaces 

We consider one dimensional examples. Multidimensional kernels can be built using tensor products of 

one- dimensional kernels. 

1. Let x G [0, 7r] and let us consider the following feature space: 

(f)(x) = (sin(a;), — p sin(2a;), —= sin(3a;), . . . , —= sin(raa;), . . .) 
V 2 v 3 v n 

Then 



00 1 1 



sin2±S 



sin^ 



K(x, y) = 4>(x) ■ 4>{x) = V] — sin(raa;) s'm(ny) = —log 

n=l 

which corresponds to the choice a n = -k=. 
2. Let x G [0,27r], h a positive number such that h < 1, and let us consider the following feature space: 

(f)(x) = (1, /&2 sin(a;), h/* cos(a;), h sin(2a;), h cos(2a;), . . . , h~ sin(raa;), h~ cos(nx), . . .) 
Then 

oo oo 

K(x, y) = 4>(x) ■ 4>{x) = 1 + 2_\ h n sin(raa;) s'm(ny) + yj h n cos(nx) cos(ny) = 

n=l n=l 

1 1-h 2 



2ir 1 — 2h cos(a; — y) + h 2 

n 

which corresponds to the choice a n = h? . 
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3. In the two examples above we have an infinite number of features, but countable. We can also construct 
cases in which the number of features is infinite and uncountable. Let us consider the following map: 



<£(x) = | ^G(s)e 8X - s | s G R c 

where G(s) is the Fourier Transform of a positive definite function, and where we work, for simplicity, 
with complex features. This corresponds to a kernel 

tf(x, y) = <£(x) • 4>(y) = f ds G(s)e^-y> s = G(x - y). 

JR d 



which corresponds to a continuum of coefficients a(s) = yG(s). 

Starting from the kernel 

Another approach consists in looking for a kernel which is known to have a representation of the form (35) 
for some set of (f>j, but whose explicit analytic form may not be known. 

In order to find a solution to this problem we need some preliminary facts. Let us call positive definite 
kernel any function ii'(x,y) on 0, X fi, with 0, C R , with the property that: 

n 

) j K(x.i,x.j)ciCj > Vxi,Xj G fi , yci,Cj G R (41) 

In the following, we will assume that 0, = [a, b] The kernel K defines an integral operator that is known 
to have a complete system of orthonormal eigenfunctions: 

dy K(x,y)0 n (y) = A n <^(x) (42) 

In 1976, Stewart [20] reported that, according to a theorem of Mercer from 1909 [13] the following state- 
ments are equivalent: 

1. The function ii'(x, y) is a positive definite kernel; 

2. 



3. The eigenvalues X n in eq. (42) are all positive; 

4. The series 



/ dxdylir(x,y)</(x)</(y)>0 Vg e C([a,b] a 

J\a,b] d 



K(*,y)= L«„i(x)^(y) 



71 = 1 



(where a 2 n = y- ) converges absolutely and uniformly. 

This leads to the following: 

Statement 2.1 Any feature vector </>(x) = (ai^>(x), ^(x), . . ., ^> n (x), . . .) such that the {a n }^ x and the 
{4>, n }'^ =1 are respectively the eigenvalues and the eigenfunctions of a positive definite kernel ii'(x, y) will 
solve problem (2.1), and the scalar product </>(x) • </>(y) has the following simple expression: 

<Kx) • <Ky) = K (^y) 

13 



A number of observations are in order: 

• Vapnik (1995) uses the condition (2) above to characterize the kernels that can be used in SVM. 
Definition (41) can be used instead, and might be more practical to work with if one has to prove the 
"admissibility" of a certain kernel. 

• There is another result similar to Mercer's one, but more general. Young (1909) proves that a kernel 
is positive definite if and only if 



/ dxdyK(x,y)g(x)g(y)>0 V</ G Xi(fi) 
Jn 



• The kernels K that can be used to represent a scalar product in the feature space are closely related 
to the theory of Reproducing Kernel Hilbert Spaces (RKHS) (see appendix A in (Girosi, 1997)[9]). In 
fact, in 1916 Moore [15] considers a more general setting for positive definite kernels, and replaces 
in eq. (41) with any abstract set E. He calls these functions positive Hermitian matrices and shows 
that to any such K one can associate a RKHS. 

In table (1) we report some commonly used kernels. 



Kernel Function 


Type of Classifier 


#(x,y) = ex p(-|l x -y|| 2 ) 


Gaussian RBF 


ii'(x,y)=(l + x-y) d 


Polynomial of degree d 


ii'(x, y) = tanh(x • y — 9) 
(only for some values of 9) 


Multi Layer Perceptron 



Table 1: Some possible kernel functions and the type of decision surface they define. 



2.4 Additional Geometrical Interpretation 

Just as Figure (2) shows why better generalization is expected from maximizing the margin, one should 

wonder: do the support vectors have any geometrical common characteristic? Are they just scattered 

points used in a linear combination? It turns out that they are not. 

In order to find the optimal decision surface, the support vector training algorithm tries to separate, as 

best as possible, the clouds defined by the data points from both classes. 

Particularly, one would expect points closer to the boundary between the classes to be more important in 

the solution than data points that are far away, since the first are harder to classify. These data points, 

in some sense, help shape and define better the decision surface than other points. Therefore, the support 

vectors are from a geometrical point of view border points. 

A direct consequence of the preceding argument delivers another important geometrical and algorithmic 

property, which is that, usually, the support vector are very few. 

These ideas can be justified algebraically through the optimality conditions derived in section 3.2.1. 

Figure (3) shows examples of the preceding geometrical interpretations with polynomial and RBF classifiers. 

2.5 An Interesting Extension: A Weighted SVM 

The original formulation of the SVM in the existing literature can be extended to handle two frequent 
cases in pattern classification and recognition: 

• An unequal proportion of data samples between the classes. 

• A need to tilt the balance or weight one class versus the other, which is very frequent when a classifi- 
cation error of one type is more expensive or undesirable than other. 
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Figure 3: Decision Surfaces given in (a) by a polynomial classifier, and in (b) by a RBF, where the 
Support Vectors are indicated in dark fill. Notice the reduce number of them and their position close to 
the boundary. In (b), the Support Vectors are the RBF centers. 

The way to derive this extension is to allow equation (37) to be: 



(43) 



where y = (yi, . . . , yi), D is a symmetric Ixl matrix with elements Dij = ytyjKfai, Xj), and C + and C~ 

are positive constants. 

Equation (18) for k = 1 now becomes: 

rain $(w,S) 



Maximize 


F(A) 


= A-l - 


-iA-M 




subject to 


A-y 

A,- 


= 
< C+l 


for y t 


= +1 




A,- 


< C~l 


for y t 


= -I 




A 


> 







— w 

2" ' 






c~ 



E e.- 



(44) 



':|/. 



and equations (19) and (20) remain unchanged. 

The quadratic program (43) can be interpreted as penalizing with higher penalty (C + or C~) the most 

undesirable type of error through equation (44) . It is also important to notice that this extension has no 

real impact on the complexity of the problem of finding the optimal vector of multipliers A, since only the 

bounding box constraints have changed. 

Notice that this extension could be changed even further to allow, for example, higher values of C for 

highly reliable or valuable data points and lower values for data points of less confidence or value. 

3 Training a Support Vector Machine 

Solving the quadratic program (37) determines the value of A* and therefore the desired decision surface 
given by equation (34). This optimization process is referred to as Training a Support Vector Machine. This 
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section covers previous, current, and possible future approaches in solving this problem. 
One important characteristic of (37) is that the quadratic form matrix D that appears in the objective 
function ( even though symmetric ) is completely dense and with size square in the number of data vectors. 
This fact implies that due to memory and computational constraints, problems with large data sets (above 
~ 5,000 samples) cannot be solved without some kind of data and problem decomposition. 
Section 3.1 deals with approaches to solving small problems, both because they constitute a natural first 
step, and also because the decomposition algorithm described in section 3.2 iteratively solves small sub- 
problems of the type given by (37). 

3.1 Previous Work 

The training of a SVM with small data sets was first approached by Vapnik et al. [5] using a constrained 
conjugate gradient algorithm. Briefly described, conjugate gradient ascent was used to explore the feasible 
region until the step would move the solution outside of it. When that happened, the largest step along 
the conjugate direction was taken, while maintaining feasibility. Every time a variable \ reached 0, 
the corresponding data point was removed (therefore reducing and approximating the solution) and the 
conjugate gradient process was re-started. 

The next approach taken by Vapnik et al. [5], was to adapt to this problem the algorithm for bounded 
large-scale quadratic programs due to More and Toraldo [16]. Originally, this algorithm uses conjugate 
gradient to explore the face of the feasible region defined by the current iterate, and gradient projection 
to move to a different face. The main modification was to only consider binding (and therefore frozen) 
those variables that were equal to one of the bounds and for which movement along the gradient would 
take them outside the feasible region. 

During the process of this research, the training problem for small data sets has been approached with 
two different algorithms and three computer packages: Zoutendijk's method of feasible directions, ( using 
CPLEX to solve the LP's ), GAMS/MINOS ( using GAMS as the modeling language and MINOS 5.4 
as the solver ), and a second-order variant of the reduced gradient method ( algorithm implemented in 
MINOS 5.4 ). A summary of these approaches and some computational results are reported next: 

3.1.1 Methods Description 

Zoutendijk's Method (case of linear constraints) [29] [1]: 

In order to solve a nonlinear problem of the form: 

Maximize /(x) 

subiect to , ,„. 

Ax <b ^ 

Ex =e 

this method follows the following skeletal approach: 

1. Find xi with A\Xi = t>i, A2X1 < D2 and Ex.\ = e, partitioning A T = [A^Aj] and b = [t>i,D2]. Let 
k = l. 

2. Given x^, A\~x.k = t>i an d A^Xk < D2, find d^, the optimal solution of: 

Maximize V/(x^) • d 
subject to 

Aid < (46) 

Ed =0 

-1 < d < 1 

3. If V/(x fc ) • d fc = 0, Stop. Else go to 4. 
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4. Find a step-size 8k, solution to: 



Maximize /(x^ + 8d) 
subject to 

< 8 < 8 max 



(47) 



where 



min i,>o(f ) if d £ ° 



d t >0\l 



00 



if d < 



with b = bo — AoXk and d = A?d 



i2Qfc- 



5. Let Xfc_|_i = Xfc + 8kdk- Let A; = A; + 1. Go to step 2. 

Step 2 involves solving a linear program, which is usually very easy to .In the case of training a SVM, step 
2 becomes: 



Maximize 


/(d) 


= (1- 


-DAk) 


•d 




subject to 


Ad 


= 










-1 < di 


< o 






for A t - = C 




< di 


< 1 






for A 8 = 




-1 < di 


< 1 






otherwise 



(48) 



and step 4 selects 6^ = min(£ p£, £ maa; ), where: 



-'opt 



d-1 -d-M 
-2dUd 



and 



mini min 

d 8 ^0 A 8 <C*,d 8 >0 



C — A 8 \ 



-A,; 



mm l — — 

A 8 >0,d,<0 di 



)} 



One interesting modification that was done to this algorithm in order to help its speed in the computer 

implementation, was to solve the problem several times with increasing upper bound C . The starting value 

of C was usually very low, and it was scaled several times until it reached the original value. The solutions 

were also scaled and used as a starting point for the following iteration. 

From a computational point of view, this method behaved a lot better than a naive constrained conjugate 

gradient implementation, both in terms of speed and graceful degradation with the increase of C. 

On the other hand, this implementation had serious difficulties in cases where most of the A 8 's were strictly 

between their bounds. The zigzagging and slow convergence it presented allowed GAMS/MINOS and 

MINOS 5.4 to outperform it by several orders of magnitude. 

GAMS/MINOS: 

GAMS is a modeling language that allows fast description and maintainability of optimization problems. 

As a language, GAMS generates the specified model and calls a user-specified solver, depending on the 

type of problem at hand. In the case of nonlinear programs, MINOS is one of these solvers. 

The work done with GAMS/MINOS was very important. At the beginning, it offered a verification of the 

implementation of Zoutendijk's method and a point of comparison in terms of speed and accuracy, but 

most important, it later pointed to the idea of using MINOS 5.4 directly, without the overhead that GAMS 

could represent. 

Another reason for considering important the work done with GAMS /MINOS was the improvement in the 

training speed due to a problem reformulation. The original problem given by (37) can be rewritten as: 
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Maximize 


F(\, 


ft) 


= A-l - 


-|A 


A, ft 










subject to 


A-y 

DA 

A 
A 




= 
= ft 

< CI 
> 





ft 



(49) 



Although strange at first sight, this transformation allows a much faster function and gradient evalua- 
tion, and was responsible for an important speedup in both steps of the solution (model generation and 
optimization). This was enough reason to use it as the formulation when using MINOS 5.4. 
MINOS 5.4: 

MINOS 5.4 solves nonlinear problems with linear constraints using Wolfe's Reduced Gradient algorithm in 
conjunction with Davidson's quasi- Newton method. Details of its implementation are described by Murtagh 
and Saunders in [17], and MINOS 5.4 User's Guide [18] and Bazaraa et al. [1] present an overview with 
some heuristics and comparisons. 

Wolfe's Reduced Gradient method depends upon reducing the dimensionality of the problem by represent- 
ing all variables in terms of an independent set of them. Under non-degeneracy assumptions (to facilitate 
this brief description), a program of the form: 

Minimize /(x) 
subject to 

Ax =b 

x > 

can be decomposed into A = [B, N] , x = (x#, xjv) with B non-singular, x# > and xjv > 0. 
By denoting the gradient V/(x) = (Vb/(x), Vjv/(x)) and a direction vector d = (d#,djv), the system 
Ad = holds for any choice of djv by letting d# = — _B _1 djv- 

Defining r = (r#, rjy) = V/(x) — Vb/(x)5 _1 A = (0, Vj\r/(x) — Vb/(x)5 _1 A) as the reduced gradient, 
it follows that V/(x) • d = rjy • djv- Therefore, in order to have a feasible direction d to be an improving 
feasible direction (feasibility and V/(x) • d < 0), a vector djv must be chosen such that rjy • djv < and 
= 0. This can be accomplished by choosing d# = — B~ 



d 3 > for 



" 1 d JV and: 



d„ 



if rj < 
if rj > 



for j G A. After determining the improving feasible direction d, a line-search procedure is used to determine 

the step-size, and an improved solution is obtained. 

Reduced gradient methods allow all components of djv to be non-zero. On the opposite side, for example, 

the simplex method for linear programming examines a similar direction-finding problem, but allow only 

one component of djv to be non-zero at a time. It is interesting to see that although the second strategy 

looks too restrictive, the first one also can result in a slow progress, as sometimes only small step-sizes are 

possible due to the fact that many components are changing simultaneously. 

In order to reach a compromise between the two strategies mentioned above, the set of non basic variables 

xjv can be further partitioned into (xs,xjv), with the corresponding decomposition of A = [S, A'] and 

d/v = (ds, djv')- The variables x s are called superbasic variables, and are intended to be the driving force 

of the iterates while xjv' is fixed and x# is adjusted to maintain feasibility [17]. 

Notice that the direction vector d can be accordingly decomposed through a linear operator Z of the form: 



" d B 




ds 


= 


d N , 





-B- X S 
I 




Zd f 



(50) 
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and now the search direction along the surface of the active constraints is characterized as being in the 
range of a matrix Z which is orthogonal to the matrix of constraint normals, i.e., 



AZ = [B,S,N'] 



-B-\S 
I 




0. (51) 



By expressing the directional derivative V/(x) • d as: 

V/(x)-d = V/(x)-Zd s = [V 5 /(x)-V B /(x)5- 1 5 , ]d 5 = r s -d s (52) 

where rs = Vs/(x) — Vb/(x)_B _1 5', and the direction finding problem can therefore be reduced to: 

Minimize rs • ds 

subject to (53) 

— ^jkil < dj < \rj\ for Xj superbasic. 

Given that the direction finding problem described by equation (53) uses a linear approximation to the 
objective function, slow and zigzagging convergence is likely to happen when the contours of / are flat or 
thin in some directions. Therefore, we would expect faster convergence when this approach is upgraded 
by taking a second-order approximation to /. More formally, the goal is to minimize a second-order 
approximation to the direction finding problem given by: 

/(x) + V/(x) • d + id • #(x)d (54) 

over the linear manifold Ad = 0. 

Using equations (52) and (50), (54) transforms into: 

min{r 5 • d 5 + id 5 • Z T H(*)Zd s } (55) 

where the matrix Z T H(x)Z is called the reduced Hessian. 

Setting the gradient of (55) equal to zero results in the system of equations: 

Z T H(x)Zd s = -r s (56) 

Once ds is available, a line-search along the direction d = Zds is performed and a new solution is obtained. 
MINOS 5.4 implements (56) with certain computational highlights [17]: 

1. During the algorithm, if it appears that no more improvement can be made with the current partition 
[B,S,N'], that is, ||rs|| < e, for a suitably chosen tolerance level e, some of the non-basic variables 
are added to the superbasics set. Using a Multiple Pricing option, MINOS allows the user to specify 
how many of them to incorporate. 

2. If at any iteration a basic or superbasic variable reaches one of its bounds, the variable is made 
non-basic. 

3. The matrices Z, -ff(x) and Z T H(x)Z are never actually computed, but are used implicitly 

4. The reduced Hessian matrix Z H{px.)Z is approximated by Br R, where R is a dense upper triangular 
matrix. 

5. A sparse LU factorization of the basis matrix B is used. 
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3.1.2 Computational Results 

In order to compare the relative speed between these methods, two different problems with small data-sets 
were solved in the same computational environment: 

1. Training a SVM with a linear classifier in the Ripley data-set. This data-set consists of 250 samples 
in two dimensions which are not linearly separable. Table 2 shows the following points of comparison: 

• The difference between GAMS/MINOS used in the original problem and in the transformed 
version (49). 

• The performance degradation suffered by the conjugate gradient implementation under the in- 
crease of the upper bound C, and on the opposite hand, the negligible effect on GAMS/MINOS 
(modified) and MINOS 5.4. 

• A considerable advantage in performance by MINOS 5.4. 

2. Training a SVM with a third degree polynomial classifier on the Sonar data-set. This data-set consists 
of 208 samples in 60 dimensions which are not linearly separable, but are polynomially separable. The 
results of this experiments are shown in Table 3 and exhibit the following points of comparison: 

• The difficulty experienced by first-order methods like Zoutendijk's method to converge when the 
values of the A 8 's are strictly between the bounds. 

• The clear advantage in solving the problem directly with MINOS 5.4, removing the overhead 
created by GAMS and incorporating the knowledge of the problem into the solution process 
through, for example, fast and exact gradient evaluation, use of symmetry in the constraint 
matrix, etc. 

• Again, a negligible effect of the upper bound C on the performance, when using MINOS. 

An important computational result is the sub-linear dependence of the training time with the dimensionality 
of the input data. In order to show this dependence, Table 4 presents the training time for randomly- 
generated 2,000 data-points problems, with different dimensionality, separability, and upper bound C . 





Methods 


c 


Conjugate Gradient 


Zoutendijk 


GAMS/MINOS 


GAMS/MINOS Modified 


MINOS 5.4 


10 


23.9 sec 


12.4 sec 


906 sec 


17.6 sec 


1.2 sec 


100 


184.1 sec 


37.9 sec 


1068 sec 


19.7 sec 


1.4 sec 


10000 


5762.2 sec 


161.5 sec 


1458 sec 


22.6 sec 


2.3 sec 



Table 2: Training time on the Ripley data-set for different methods and upper bound C. GAMS/MINOS 
Modified corresponds to the reformulated version of the problem. 





Methods 


c 


Zoutendijk 


GAMS/MINOS Modified 


MINOS 5.4 


10 


4381.2 sec 


67.0 sec 


3.3 sec 


100 


N/A 


67.1 sec 


3.3 sec 


10000 


N/A 


67.1 sec 


3.3 sec 



Table 3: Training time on the Sonar Dataset for different methods and upper bound C. 



3.2 A New Approach to Large Database Training 

As mentioned before, training a SVM using large data sets (above ~ 5,000 samples) is a very difficult 
problem to approach without some kind of data or problem decomposition. To give an idea of some 
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memory requirements, an application like the one described later in section 3 involves 50,000 training 

samples, and this amounts to a quadratic form whose matrix D has 2.5 • 10 9 entries that would need, using 

an 8-byte floating point representation, 20,000 Megabytes = 20 Gigabytes of memory! 

In order to solve the training problem efficiently, we take explicit advantage of the geometric interpretation 

introduced in Section 2.4, in particular, the expectation that the number of support vectors will be very 

few. If we consider the quadratic programming problem given by (37), this expectation translates into 

having many of the components of A equal to zero. 

In order to decompose the original problem, one can think of solving iteratively the system given by (37), 

but keeping fixed at zero level, those components \ associated with data points that are not support 

vectors, and therefore only optimizing over a reduced set of variables. 

To convert the previous description into an algorithm we need to specify: 

1. Optimality Conditions: These conditions allow us to decide computationally, if the problem has 
been solved optimally at a particular global iteration of the original problem. Section 3.2.1 states and 
proves optimality conditions for the QP given by (37). 

2. Strategy for Improvement: If a particular solution is not optimal, this strategy defines a way to 
improve the cost function and is frequently associated with variables that violate optimality conditions. 
This strategy will be stated in section 3.2.2. 

After presenting optimality conditions and a strategy for improving the cost function, section 3.2.3 intro- 
duces a decomposition algorithm that can be used to solve large database training problems, and section 
3.2.4 reports some computational results obtained with its implementation. 

3.2.1 Optimality Conditions 

In order to be consistent with common standard notation for nonlinear optimization problems, the 
quadratic program (37) can be rewritten in minimization form as: 



Minimize 

A 


W(A) 


= -A-1 + ^A-DA 


subject to 


A-y 
A-Cl 


= (/x) 
< (Y 




-A 


<o (n 



(57) 



where /i, Y = (v\, . . . , vi) and II = (tti, . . . , iri) are the associated Kuhn- Tucker multipliers. 

Since D is a positive semi-definite matrix (see end of section 2.3.3) and the constraints in (57) are linear, 

the Kuhn- Tucker (KT) conditions are necessary and sufficient for optimality, and they are: 





Dimension 




Separable 


Non- Separable 


c 


4 


16 


256 


4 


16 


256 


10 


60.7 sec 


106.4 sec 


613.5 sec 


292.9 sec 


476.0 sec 


1398.2 sec 


100 


36.0 sec 


69.2 sec 


613.7 sec 


313.5 sec 


541.0 sec 


2369.4 sec 


10000 


21.8 sec 


56.2 sec 


623.0 sec 


327.4 sec 


620.6 sec 


3764.1 sec 



Table 4: Training time on a Randomly-generated Dataset for different dimensionality and upper bound C. 
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+ Y- 


n + 


^y 


= 


Y-(A- 


-Cl) 
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n a 








= 


Y 








> o 


n 








> o 


A-y 








= 


A -CI 








< o 


-A 
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In order to derive further algebraic expressions from the optimality conditions(58) , we assume the existence 
of some A,- such that < A 8 - < C (see end of section 2.3.2), and consider the three possible values that each 
component of A can have: 

1. Case: < A; < C: 

From the first three equations of the KT conditions we have: 

(DA)i -l + fi yi = (59) 

Noticing that 

i i 

(DA), = Y^ ^jyjViK(y.i, xj) = yi ^ Aj2/jir(xi, xj) 

3=1 3=1 

and that for < A 8 - < C, 

i i 

/(xi) = sign(^2\jyjK(x.i,x.j) + b) = ^ XjyjKfa,^) + b = y % 

3=1 3=1 

we obtain the following: 

(DA),- = ^ = 1 - - (60) 

Vi Vi 

By substituting (60) into (59) we finally obtain that 

fi = b (61) 

Therefore, at an optimal solution A*, the value of the multiplier fj, is equal to the optimal threshold 
b*. 

2. Case: A 8 = C: 

From the first three equations of the KT conditions we have: 

(DA) t -l + vi + nyi = (62) 

By defining 

£ 

S(*i) = Yl A J%'-^( x i' X j) + b 
3=1 
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and noticing that 



i 
{DA)i = yi Y^ ^jVjK^i, xj) = yigfa) - y x b 



equation (62) can be written as: 



Vig{xi) - y t b -1 + v t + fiyi = 



By combining fj, = b (derived from case 1) and requiring V{ > we finally obtain: 

Vidi^i) < 1 (63) 

3. Case: A 8 - = 0: 

From the first three equations of the KT conditions we have: 

(M), - 1 - 7T t - + HVi = (64) 

By applying a similar algebraic manipulation as the one described for case 2, we obtain 

Vidi^i) > 1 (65) 

3.2.2 Strategy for Improvement 

In order to incorporate the optimality conditions and the expectation that most A 8 's will be zero into an 
algorithm, we need to derive a way to improve the objective function value using this information. To do 
this, let us decompose A in two vectors A# and Ajv, where Ajv = 0, and B and N partition the index 
set, and that the optimality conditions hold in the subproblem defined only for the variables in B. In 
further sections, the set B will be referred to as the working set. Under this decomposition the following 
statements are clearly true: 

• We can replace \ = 0, i G B, with Xj = 0, j G N, without changing the cost function or the feasibility 
of both the subproblem and the original problem. 

• After such a replacement, the new subproblem is optimal if and only if yjg(xj) > 1. This follows 
from equation (65) and the assumption that the subproblem was optimal before the replacement was 
done. 

The previous statements suggest that replacing variables at zero levels in the subproblem, with variables 
Xj = 0, j G N that violate the optimality condition yjg(xj) > 1, yields a subproblem that, when 
optimized, improves the cost function while maintaining feasibility. The following proposition states this 
idea formally. 

Proposition: Given an optimal solution of a subproblem defined on B, the operation of replacing A 8 - = 0, 
i G B, with Xj = 0, j G A, satisfying yjg(xj) < 1 generates a new subproblem that when optimized, 
yields a strict improvement of the objective function W(A). 

Proof: We assume again the existence of X p such that < X p < C. Let us also assume without loss of 
generality that y p = yj (the proof is analogous if y p = —yj). Then, there is some e > such that X p — S > 0, 
for 8 G (0, e). Notice also that #(x p ) = y p . Now, consider A = A + 8ej — 8e p , where ej and e p are the jth 
and pth unit vectors, and notice that the pivot operation can be handled implicitly by letting 8 > and 
by holding A 8 - = 0. The new cost function VT(A) can be written as: 
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W(A) = -A-1 + -A-M 



-A • 1 + - [A • DA + 2A • D(8e 3 - 8e p ) + (8e 3 - 8e p ) ■ D(8e 3 - 8e p ) 



W(A) + 8 



g( x j) ~ b 1 ! b 
8 2 



s 2 

— [R (xj, xj) + R (x p , x p ) - 2y p yjR (x p , Xj) 



W(A) + 8 [g(x 3 )y 3 - 1] + — [R(x 3 ,x 3 ) + R(x p ,x p ) - 2y p y 3 R(x p ,x 3 )] 

Therefore, since g(xj)yj < 1, by choosing 8 small enough we have T'F(A) < W(A). 
q.e.d 

3.2.3 The Decomposition Algorithm 

Suppose we can define a fixed-size working set B, such that \B\ < £, and it is big enough to contain all 
support vectors (A 8 - > 0), but small enough such that the computer can handle it and optimize it using 
some solver. Then the decomposition algorithm can be stated as follows: 

1. Arbitrarily choose \B\ points from the data set. 

2. Solve the subproblem defined by the variables in B. 

3. While there exists some j £ N , such that g(xj)yj < 1, where 

i 

P =i 

replace A 8 - = 0, i G B, with Xj = and solve the new subproblem. 

Notice that this algorithm will strictly improve the objective function at each iteration and therefore will 
not cycle. Since the objective function is bounded (W(A) is convex and quadratic, and the feasible region 
is bounded), the algorithm must converge to the global optimal solution in a finite number of iterations. 
Figure 4 gives a geometric interpretation of the way the decomposition algorithm allows the redefinition 
of the separating surface by adding points that violate the optimality conditions. 

3.2.4 Computational Implementation and Results 

We have implemented the decomposition algorithm using the transformed problem defined by equation 
(49) and MINOS 5.4 as the solver. 

Notice that the decomposition algorithm is rather flexible about the pivoting strategy, that is, the way it 
decides which and how many new points to incorporate into the working set B. Our implementation uses 
two parameters to define the desired strategy: 

• Lookahead: this parameter specifies the maximum number of data points the pricing subroutine 
should use to evaluate optimality conditions (Case 3). If Lookahead data points have been examined 
without finding a violating one, the subroutine continues until it finds the first one , or until all data 
points have been examined. In the latter case, global optimality has been obtained. 

• Newlimit: this parameter limits the number of new points to be incorporated into the working set 
B. 
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Figure 4: (a) A sub-optimal solution where the non-filled points have A = but are violating optimality 
conditions by being inside the ±1 area, (b) The decision surface gets redefined. Since no points with A = 
are inside the ±1 area, the solution is optimal. Notice that the size of the margin has decreased, and the 
shape of the decision surface has changed. 

The computational results that we present in this section have been obtained using real data from our Face 
Detection System, which is described in Section 4. 

Figure 5 shows the training time and the number of support vectors obtained when training the system 
with 5,000, 10,000, 20,000, 30,000, 40,000, 49,000, and 50,000 data points. We must emphasize that the 
last 1,000 data points were collected in the last phase of bootstrapping of the Face Detection System, and 
therefore make the training process harder, since they correspond to errors obtained with a system that 
was already very accurate. Figure 6 shows the relationship between the training time and the number 
of support vectors, as well as the number of global iterations (the number of times the decomposition 
algorithm calls the solver). Notice the smooth relation between the number of support vectors and the 
training time, and the jump from 11 to 15 global iterations as we go from 49,000 to 50,000 samples. This 
increase is responsible for the increase in the training time. The system, using a working set of 1200 
variables was able to solve the 50,000 data points problem using only 25Mb of RAM. 

Figure 7 shows the effect on the training time due to the parameter Newlimit and the size of the working 
set, using 20,000 data points. Notice the clear improvement as Newlimit is increased. This improvement 
suggests that in some way, the faster new violating data points are brought into the working set, the faster 
the decision surface is defined, and optimality is reached. Notice also that, if the working set is too small or 
too big compared to the number of support vectors (746 in the case of 20,000 samples), the training time 
increases. In the first case, this happens because the algorithm can only incorporate new points slowly, 
and in the second case, this happens because the solver takes longer to solve the subproblem as the size of 
the working set increases. 

3.3 Improving the Training of SVM: Future Directions 

The algorithm described in Section 3.2.3 suggests two main areas where improvements can be made through 
future research. These two areas are: 
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Figure 5: (a) Training Time on a SPARCstation-20. (b) Number of Support Vectors obtained after Training 

1. The Solver: The second-order variant of the reduced gradient method implemented by MINOS 5.4 
has given very good results so far in terms of accuracy, robustness and performance. However, this 
method is a general nonlinear optimization method that is not designed in particular for quadratic 
programs, and in the case of SVM's, is not designed in particular for the special characteristics of the 
problem. Having as a reference the experience obtained with MINOS 5.4, new approaches to a tailored 
solver through, for example, projected Newton [2] or interior point methods [7], should be attempted. 

At this point it is not clear whether the same type of algorithm is appropriate for all stages of the 
solution process. To be more specific, it could happen that an algorithm performs well with few 
non-zero variables at early stages, and then is outperformed by others when the number of non-zero 
variables reaches some threshold. In particular, we learned that the number of non-zero variables that 
satisfy < A 8 - < C has an important effect on the performance of the solver. 

2. The Pivoting Strategy: This area offers great potential for improvements through some ideas we 
plan to implement. The improvements are based on some qualitative characteristics of the training 
process that have been observed: 

• During the execution of the algorithm, as much as 40% of the computational effort is dedicated 
to the evaluation of the optimality conditions. At final stages, it is common to have all the data 
points evaluated, yet only to collect very few of them to incorporate them to the working set. 

• Only a small portion of the input data is ever brought into the working set (about 16% in the 
case of the face detection application). 

• Out of the samples that ever go into the working set, about 30% of them enter and exit this set 
at least once. These vectors are responsible for the first characteristic mentioned above. 

Possible future strategies that exploit these characteristics are: 

• Keep a list or file with all or part of the input vectors that have exited the working set. At 
the pricing stage, when the algorithm computes the optimality conditions, evaluate these data 
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Figure 6: (a) Number of Support Vectors versus Training Time on a SPARCstation-20. Notice how the 
Number of support vectors is a better indicator of the increase in training time than the number of samples 
alone, (b) Number of global iterations performed by the algorithm. Notice the increase experimented when 
going from 49000 to 50000 samples. This increase in the number of iterations is responsible for the increase 
in the training time 

poits before other data points to determine the entering vectors. This strategy is analogous 
to one sometimes used in the revised simplex method where the algorithm keeps track of the 
basic variables that have become non-basic. In the case of training of SVM's, the geometric 
interpretation of this heuristic is to think that if a point was a support vector at some iteration, 
it was more or less close to the boundary between the classes, and as this boundary is refined or 
fine-tuned, it is possible for it to switch from active to non-active several times. This heuristic 
could be combined with main memory and cache management policies used in computer systems. 

• During the pricing stage, instead of bringing into the working set the first k points that violate 
optimality conditions, we could try to determine r violating data points, with r > k and choose 
from these the k most violating points. This is done under the geometric idea that the most 
violating points help in defining the decision surface faster and therefore save time in future 
iterations. 

• These last two approaches can be combined by keeping track not only of the points exiting the 
working set, but also of the remaining violating data points as well. 

So far in the description and implementation of the decomposition algorithm, we have assumed that enough 
memory is available to solve a working set problem that contains all of the support vectors. However, 
some applications may require more support vectors than the available memory can manage. One possible 
approach that can be taken is to approximate the optimal solution by the best solution that can be obtained 
with the current working set size. The present algorithm and implementation can be easily extended to 
handle this situation by replacing support vectors with < A 8 - < C with new data points. Other more 
complex approaches that can be pursued to obtain an optimal solution should be the subject of future 
research. 
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Figure 7: (a) Training time for 20,000 samples with different values of Newlimit, using a working set of 
size fOOO and Lookahead=f 0,000. (b) Training time for 20,000 samples with different sizes of the working 
set, using Newlimit=size of the working set, and Lookahead=f 0,000. 

4 SVM Application: Face Detection in Images 

This section introduces a Support Vector Machine application for detecting vertically oriented and unoc- 

cluded frontal views of human faces in grey level images. It handles faces over a wide range of scales and 

works under different lighting conditions, even with moderately strong shadows. 

The face detection problem can be defined as follows: Given as input an arbitrary image, which could be 

a digitized video signal or a scanned photograph, determine whether or not there are any human faces in 

the image, and if there are, return an encoding of their location. The encoding in this system is to fit each 

face in a bounding box defined by the image coordinates of the corners. 

Face detection as a computer vision task has many applications. It has direct relevance to the face 

recognition problem, because the first important step of a fully automatic human face recognizer is usually 

identifying and locating faces in an unknown image. Face detection also has potential application in 

human-computer interfaces, surveillance systems, census systems, etc. 

From the standpoint of this paper, face detection is interesting because it is an example of a natural and 

challenging problem for demonstrating and testing the potentials of Support Vector Machines. There are 

many other object classes and phenomena in the real world that share similar characteristics, for example, 

tumor anomalies in MRI scans, structural defects in manufactured parts, etc. A successful and general 

methodology for finding faces using SVM's should generalize well for other spatially well-defined pattern 

and feature detection problems. 

It is important to remark that face detection, like most object detection problems, is a difficult task due 

to the significant pattern variations that are hard to parameterize analytically. Some common sources of 

pattern variations are facial appearance, expression, presence or absence of common structural features, 

like glasses or a moustache, light source distribution, shadows, etc. 

This system works by testing candidate image locations for local patterns that appear like faces using 

a classification procedure that determines whether or not a given local image pattern is a face or not. 
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Therefore, the face detection problem is approached as a classification problem given by examples of 2 
classes: faces and non-faces. 

4.1 Previous Systems 

The problem of face detection has been approached with different techniques in the last few years. This 

techniques include Neural Networks [4] [19], detection of face features and use of geometrical constraints 

[27], density estimation of the training data [14], labeled graphs [12] and clustering and distribution-based 

modeling [22] [21]. 

Out of all these previous works, the results of Sung and Poggio [22][21], and Rowley et al. [19] reflect 

systems with very high detection rates and low false positive detection rates. 

Sung and Poggio use clustering and distance metrics to model the distribution of the face and non-face 

manifold, and a Neural Network to classify a new pattern given the measurements. The key of the quality 

of their result is the clustering and use of combined Mahalanobis and Euclidean metrics to measure the 

distance from a new pattern and the clusters. Other important features of their approach is the use of 

non-face clusters, and the use of a bootstrapping technique to collect important non-face patterns. One 

drawback of this technique is that it does not provide a principled way to choose some important free 

parameters like the number of clusters it uses. 

Similarly, Rowley et al. have used problem information in the design of a retinally connected Neural 

Network that is trained to classify faces and non-faces patterns. Their approach relies on training several 

NN emphasizing subsets of the training data, in order to obtain different sets of weights. Then, different 

schemes of arbitration between them are used in order to reach a final answer. 

The approach to the face detection system with a SVM uses no prior information in order to obtain the 

decision surface, this being an interesting property that can be exploited in using the same approach for 

detecting other objects in digital images. 

4.2 The SVM Face Detection System 

This system, as it was described before, detects faces by exhaustively scanning an image for face-like 

patterns at many possible scales, by dividing the original image into overlapping sub-images and classifying 

them using a SVM to determine the appropriate class, that is, face or non-face. Multiple scales are handled 

by examining windows taken from scaled versions of the original image. 

Clearly, the major use of SVM's is in the classification step, and it constitutes the most critical and 

important part of this work. Figure 9 gives a geometrical interpretation of the way SVM's work in the 

context of face detection. 

More specifically, this system works as follows: 

1. A database of face and non-face 19x19 pixel patterns, assigned to classes +1 and -1 respectively, is 
trained on, using the support vector algorithm. A 2nd-degree polynomial kernel function and an upper 
bound C = 200 are used in this process obtaining a perfect training error. 

2. In order to compensate for certain sources of image variation, some preprocessing of the data is 
performed: 

• Masking: A binary pixel mask is used to remove some pixels close to the boundary of the window 
pattern allowing a reduction in the dimensionality of the input space from 19x19=361 to 283. 
This step is important in the reduction of background patterns that introduce unnecessary noise 
in the training process. 

• Illumination gradient correction: A best-fit brightness plane is subtracted from the unmasked 
window pixel values, allowing reduction of light and heavy shadows. 

• Histogram equalization: A-histogram equalization is performed over the patterns in order to 
compensate for differences in illumination brightness, different cameras response curves, etc. 
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3. Once a decision surface has been obtained through training, the run-time system is used over images 
that do not contain faces, and misclassifications are stored so they can be used as negative examples 
in subsequent training phases. Images of landscapes, trees, buildings, rocks, etc., are good sources 
of false positives due to the many different textured patterns they contain. This bootstrapping step, 
which was successfully used by Sung and Poggio [22] is very important in the context of a face detector 
that learns from examples because: 

• Although negative examples are abundant, negative examples that are useful from a learning 
point of view are very difficult to characterize and define. 

• By approaching the problem of object detection, and in this case of face detection, by using the 
paradigm of binary pattern classification, the two classes, object and non-object are not equally 
complex since the non-object class is broader and richer, and therefore needs more examples in 
order to get an accurate definition that separates it from the object class. Figure 8 shows an image 
used for bootstrapping with some misclassifications, that were later used as negative examples. 

4. After training the SVM, we incorporate it as the classifier in a run-time system very similar to the 
one used by Sung and Poggio [22] [21] that performs the following operations: 

• Re-scale the input image several times. 

• Cut 19x19 window patterns out of the scaled image. 

• Preprocess the window using masking, light correction and histogran equalization. 

• Classify the pattern using the SVM. 

• If the class corresponds to a face, draw a regtangle aroung the face in the output image. 



4.2.1 Experimental Results 

To test the run-time system, we used two sets of images. The set A, contained 313 high-quality images 
with same number of faces. The set B, contained 23 images of mixed quality, with a total of 155 faces. 
Both sets were tested using our system and the one by Sung and Poggio [22] [21]. In order to give true 
meaning to the number of false positives obtained, it is important to state that set A involved 4,669,960 
pattern windows, while set B 5,383,682. Table 5 shows a comparison between the 2 systems. 
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Detection Rate 


False Detections 


Detection Rate 


False Detections 


Ideal System 


100% 





100% 





SVM 


97.12 % 


4 


74.19% 


20 


Sung and Poggio 


94.57 % 


2 


74.19% 


11 



Table 5: Performance of the SVM face detection system 

Figures 10, 11, 12, 13, 14 and 15 present some output images of our system. These images were not used 
during the training phase of the system. 

4.3 Future Directions in Face Detection and SVM Applications 

Future research in SVM application can be divided into three main categories or topics: 

1. Simplification of the SVM: One drawback for using SVM in some real-life applications is the large 
number of arithmetic operations that are necessary to classify a new input vector. Usually, this number 
is proportional to the dimension of the input vector and the number of support vectors obtained. In 
the case of face detection, for example, this is ~ 283,000 multiplications per pattern! 

The reason behind this overhead is in the roots of the technique, since SVM's define the decision 
surface by explicitly using data points. This situation causes a lot of redundancy in most cases, and 
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can be solved by relaxing the constraint of using data points to define the decision surface. This is a 
topic of current research conducted by Burges [6] at Bell Laboratories, and it is of great interest for us, 
because in order to simplify the set of support vectors, one needs to solve highly-nonlinear constrained 
optimization problems. 

Since a closed form solution exists for the case of kernel functions that are 2nd. degree polynomials, 
we are using a simplified SVM [6] in our current experimental face detection system that gains an 
acceleration factor of 20, without degrading the quality of the classifications. 

2. Detection of other objects: We are interested in using SVM's to detect other objects in digital 
images, like cars, airplanes, pedestrians, etc. Notice that most of these objects have very different 
appearance, depending on the viewpoint. In the context of face detection, an interesting extension 
that could lead to better understanding and approach to other problems, is the detection of tilted 
and rotated faces. It is not clear at this point, whether these different view of the same object can be 
treated with a single classifier, or if they should be treated separately. 

3. Use of multiple classifiers: The use of multiple classifiers offers possibilities that can be faster 
and/or more accurate. Rowley et al. [19] have successfully combined the output from different neural 
networks by means of different schemes of arbitration in the face detection problem. Sung and Poggio 
[22] [21] use a first classifier that is very fast as a way to quickly discard patterns that are clearly 
non-faces. These two references are just examples of the combination of different classifiers to produce 
better systems. 

Our current experimental face detection system performs an initial quick-discarding step using a SVM 
trained to separate clearly non-faces from probable faces using just 14 averages taken from different 
areas of the window pattern. This classification can be done about 300 times faster and is currently 
discarding more than 99% of input patterns. More work will be done in the near future in this area. 

The classifiers to be combined do not have to be of the same kind. An interesting type of classifier 
that we will consider is Discriminant Adaptative Nearest Neighbor due to Hastie et al. [11] [10]. 

5 Conclusions 

In this paper we presented a novel decomposition algorithm to train Support Vector Machines. We have 
successfully implemented this algorithm and solved large dataset problems using acceptable amounts of 
computer time and memory. Work in this area can be continued, and we are currently studying new 
techniques to improve the performance and quality of the training process. 

The Support Vector Machine is a very new technique, and, as far as we know, this paper presents the 
second problem- solving application to use SVM, after Vapnik et al. [5, 8, 24] used it in the character 
recognition problem, in 1995. 

Our Face Detection System performs as well as other state-of-the-art systems, and has opened many 
interesting questions and possible future extensions. From the object detection point of view, our ultimate 
goal is to develop a general methodology that extends the results obtained with faces to handle other 
objects. From a broader point of view, we also consider interesting the use of the function approximation 
or regression extension that Vapnik [24] has done with SVM, in many different areas where Neural Networks 
are currently used. 

Finally, another important contribution of this paper is the application of OR-based techniques to domains 
like Statistical Learning and Artificial Intelligence. We believe that in the future, other tools like duality 
theory, interior point methods and other optimization techniques and concepts, will be useful in obtaining 
better algorithms and implementations with a solid mathematical background. 
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Figure 8: Some false detections obtained with the first version of the system. This false positives were 
later used as negative examples ( class -1 ) in the training process 
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Figure 9: Geometrical Interpretation of how the SVM separates the face and non-face classes. The patterns 
are real support vectors obtained after training the system. Notice the small number of total support vectors 
and the fact that a higher proportion of them correspond to non-faces. 
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Figure 10: Faces 
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Figure 11: Faces 
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Figure 12: Faces 
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Figure 15: Faces 
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